# This script analyzes the results of LC_processVS_hpc.R
Packages <- c("rstan", "coda", "ggmcmc", "bayesplot", "tidyverse",
"loo", "purrr")
suppressMessages(invisible(lapply(Packages, library, character.only=TRUE)))
theme_set(theme_bw())
setwd("C:/Users/tms1044/Desktop/lc_mod")
n.cores <- 4
chn.dir <- "out/"
sum.dir <- "summaries/"
mods <- c("Cens", "Clim", "ClimCens", "Topo",
"TopoCens", "TopoClim", "TopoClimCens")
# GGMCMC
ac_cols <- c(rgb(0,0,0,0.5), rgb(1,0,0,0.5), rgb(0,1,0,0.5),
rgb(0,0,1,0.5), rgb(1,1,0,0.5), rgb(1,0,1,0.5))
beta.fls <- list.files(sum.dir, pattern="out_betas_", full.names=TRUE)
names(beta.fls) <- mods
map(beta.fls, readRDS) %>%
map(., ~(ggs_density(ggs(.x)) + facet_wrap(~Parameter, scales="free")))
## $Cens
##
## $Clim
##
## $ClimCens
##
## $Topo
##
## $TopoCens
##
## $TopoClim
##
## $TopoClimCens
map(beta.fls, readRDS) %>%
map(., ~(ggs_traceplot(ggs(.x)) + facet_wrap(~Parameter, scales="free")))
## $Cens
##
## $Clim
##
## $ClimCens
##
## $Topo
##
## $TopoCens
##
## $TopoClim
##
## $TopoClimCens
map(beta.fls, readRDS) %>% map(., ~(ggs_autocorrelation(ggs(.x)) +
facet_wrap(~Parameter) +
scale_fill_manual(values=rep(NA,6)) +
scale_colour_manual(values=ac_cols)))
## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
## $Cens
##
## $Clim
##
## $ClimCens
##
## $Topo
##
## $TopoCens
##
## $TopoClim
##
## $TopoClimCens
map(beta.fls, readRDS) %>% map(., ~(ggs_geweke(ggs(.x))))
## $Cens
##
## $Clim
##
## $ClimCens
##
## $Topo
##
## $TopoCens
##
## $TopoClim
##
## $TopoClimCens
map(beta.fls, readRDS) %>% map(., ~(ggs_Rhat(ggs(.x))))
## $Cens
##
## $Clim
##
## $ClimCens
##
## $Topo
##
## $TopoCens
##
## $TopoClim
##
## $TopoClimCens
map(beta.fls, readRDS) %>%
map(., ~(ggs_pairs(ggs(.x), lower=list(continuous="density"))))
## $Cens
##
## $Clim
##
## $ClimCens
##
## $Topo
##
## $TopoCens
##
## $TopoClim
##
## $TopoClimCens
map(beta.fls, readRDS) %>%
map(., ~(ggs_running(ggs(.x)) + facet_wrap(~Parameter, scales="free")))
## $Cens
##
## $Clim
##
## $ClimCens
##
## $Topo
##
## $TopoCens
##
## $TopoClim
##
## $TopoClimCens
map(beta.fls, readRDS) %>%
map(., ~(ggs_compare_partial(ggs(.x)) + facet_wrap(~Parameter, scales="free")))
## $Cens
##
## $Clim
##
## $ClimCens
##
## $Topo
##
## $TopoCens
##
## $TopoClim
##
## $TopoClimCens
# effective sample size
par(mfrow=c(3,3))
map(beta.fls, readRDS) %>%
walk(., ~hist(effectiveSize(.), main="", xlab="Effective Sample Size"))
par(mfrow=c(3,3))
map(beta.fls, readRDS) %>%
walk(., ~hist(effectiveSize(.)/(250*6), main="",
xlab="n_eff/n_iter", xlim=c(0,1.5)))
par(mfrow=c(1,1))
# WAIC
waic.ls <- map(list.files(sum.dir, pattern="waic", full.names=TRUE), readRDS)
names(waic.ls) <- mods
map_dbl(waic.ls, ~(.$waic)) %>% sort
## Topo TopoCens TopoClimCens TopoClim Clim
## -796749.0 -795971.0 -785065.4 -772394.8 -767220.5
## ClimCens Cens
## -763440.7 -759462.5
compare(waic.ls[[1]], waic.ls[[2]], waic.ls[[3]], waic.ls[[4]],
waic.ls[[5]], waic.ls[[6]], waic.ls[[7]])
## waic se_waic elpd_waic se_elpd_waic p_waic
## waic.ls[[4]] -796749.0 1075.3 398374.5 537.6 84739.9
## waic.ls[[5]] -795971.0 1105.9 397985.5 552.9 84307.0
## waic.ls[[7]] -785065.4 1116.2 392532.7 558.1 77443.9
## waic.ls[[6]] -772394.8 1087.3 386197.4 543.7 80668.8
## waic.ls[[2]] -767220.5 1082.6 383610.2 541.3 83730.8
## waic.ls[[3]] -763440.7 1110.2 381720.4 555.1 84788.8
## waic.ls[[1]] -759462.5 1104.8 379731.2 552.4 81700.1
## se_p_waic
## waic.ls[[4]] 24.9
## waic.ls[[5]] 28.6
## waic.ls[[7]] 26.5
## waic.ls[[6]] 24.3
## waic.ls[[2]] 24.9
## waic.ls[[3]] 27.3
## waic.ls[[1]] 29.6
# LOO
loo.ls <- map(list.files(sum.dir, pattern="loo", full.names=TRUE), readRDS)
names(loo.ls) <- mods
map_dbl(loo.ls, ~(.$looic)) %>% sort
## Topo TopoCens TopoClimCens TopoClim Clim
## -753354.2 -751749.9 -740253.9 -728311.2 -723893.4
## ClimCens Cens
## -720686.0 -715829.4
compare(loo.ls[[1]], loo.ls[[2]], loo.ls[[3]], loo.ls[[4]],
loo.ls[[5]], loo.ls[[6]], loo.ls[[7]])
## looic se_looic elpd_loo se_elpd_loo p_loo se_p_loo
## loo.ls[[4]] -753354.2 1080.4 376677.1 540.2 106437.3 57.2
## loo.ls[[5]] -751749.9 1110.6 375875.0 555.3 106417.5 58.0
## loo.ls[[7]] -740253.9 1120.6 370126.9 560.3 99849.7 57.4
## loo.ls[[6]] -728311.2 1091.6 364155.6 545.8 102710.6 56.6
## loo.ls[[2]] -723893.4 1086.5 361946.7 543.3 105394.3 57.5
## loo.ls[[3]] -720686.0 1115.1 360343.0 557.5 106166.2 58.1
## loo.ls[[1]] -715829.4 1108.8 357914.7 554.4 103516.6 58.1
map(loo.ls, plot)
## $Cens
## NULL
##
## $Clim
## NULL
##
## $ClimCens
## NULL
##
## $Topo
## NULL
##
## $TopoCens
## NULL
##
## $TopoClim
## NULL
##
## $TopoClimCens
## NULL
# OOS prediction
mspe.ls <- map(list.files(sum.dir, pattern="MSPE", full.names=TRUE), readRDS)
names(mspe.ls) <- mods
map_dbl(mspe.ls, sqrt) %>% sort
## Clim ClimCens Cens TopoClim TopoClimCens
## 0.1810194 0.1811403 0.1854589 0.1893665 0.1899082
## Topo TopoCens
## 0.1909966 0.1914862
lpd.ls <- map(list.files(sum.dir, pattern="lpd", full.names=TRUE), readRDS)
names(lpd.ls) <- mods
map_dbl(lpd.ls, c) %>% sort
## Clim ClimCens Topo Cens TopoCens
## -2.170384 -2.139518 -2.103142 -2.101887 -2.082612
## TopoClim TopoClimCens
## -2.031828 -2.006136
# Gelman
par(mfrow=c(3,3))
map(list.files(sum.dir, pattern="gelman", full.names=TRUE), readRDS) %>%
walk2(.x=., .y=mods, ~({hist(x=.$psrf[,2], main=.y, xlim=c(0.95, 1.2));
abline(v=1.1, lty=3, col="red")}))
map(list.files(sum.dir, pattern="out_betas", full.names=TRUE), readRDS) %>%
walk2(.x=., .y=mods, ~(gelman.plot(x=.x, ylab=paste(.y, "shrink factor"),
ylim=c(0.95, 1.2))))
# Geweke
gew.ls <- map(list.files(sum.dir, pattern="geweke", full.names=TRUE), readRDS)
for(m in 1:length(mods)) {
par(mfrow=c(2,3))
for(i in 1:length(gew.ls[[m]])) {
plot(density(gew.ls[[m]][[i]]$z), main=mods[m], xlim=c(-5,5), ylim=c(0,0.5))
curve(dnorm(x, 0, 1), from=-5, to=5, add=TRUE, col="dodgerblue")
abline(v=c(-1.96, 1.96), lty=3)
}
}
map(list.files(sum.dir, pattern="out_betas", full.names=TRUE),
readRDS) %>% walk(geweke.plot)
